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Abstract. We study the phase diagram of the standard pair approximation equations for two different 
models in population dynamics, the susceptible-infective-recovered-susceptible model of infection spread 
and a predator-prey interaction model, on a network of homogeneous degree k. These models have similar 
phase diagrams and represent two classes of systems for which noisy oscillations, still largely unexplained, 
are observed in nature. We show that for a certain range of the parameter k both models exhibit an 
oscillatory phase in a region of parameter space that corresponds to weak driving. This oscillatory phase, 
however, disappears when k is large. For k = 3,4, we compare the phase diagram of the standard pair 
approximation equations of both models with the results of simulations on regular random graphs of the 
same degree. We show that for parameter values in the oscillatory phase, and even for large system sizes, 
the simulations either die out or exhibit damped oscillations, depending on the initial conditions. We 
discuss this failure of the standard pair approximation model to capture even the qualitative behavior 
of the simulations on large regular random graphs and the relevance of the oscillatory phase in the pair 
approximation diagrams to explain the cycling behavior found in real populations. 



PACS. 87.23.Cc Population dynamics and ecological pattern formation 
87.10.Rt Monte Carlo simulations 



87.10.Ca Analytical theories 



1 Introduction 

A common paradigm in population dynamics studies is 
to assume that population is not spatially distributed so 
that individuals mix perfectly and contact each other with 
equal probability. Thus, in the limit of infinite population 
the time evolution of the system is described in terms of 
, the densities of different subpopulations and governed by 
a set of differential equations which can be deduced from 
\ the law of mass action |1I2) . Another approach is to use 
stochastic dynamics on lattices or on more general graphs 
' where the discrete variables associated to each node rep- 
resent the presence and attributes of an individual, e.g. 
susceptible, infective or recovered state in the case of epi- 
demic dynamics |3I4) . These studies have shown that spa- 
tial correlations neglected in the perfect mixing assump- 
tion play an important role in the behavior of population 
dynamics on graphs, and therefore also in real popula- 
tions. 

In order to try to capture the effects of spatial cor- 
relations in the framework of an analytic description the 
standard pair approximation (PA) as well as several im- 
provements based on including higher order correlations or 
on different closure assumptions at the level of pairs have 
been proposed in the context of ecological and epidemio- 
logical models |ni6l7l8l9ll()lllll2ll.Sll4HKIlriTT7] . The per- 



formance of the PA was compared with the perfect mixing 
or mean field approximation (MFA) for several different 
population dynamics models on lattices and other graphs 
[14115116117118] . In these studies, the PA was found to per- 
form much better than the MFA in the prediction of the 
steady state of the system, and to give a good quantita- 
tive agreement in a large region of parameter space. On 
the other hand, a recent study of an epidemiological model 
has shown that an overall accurate prediction may require 
going beyond the PA to higher order cluster approxima- 
tions [TO]. 

While being more accurate than the MFA, the PA 
equations are still simple enough to allow for a complete 
description of their phase diagram, which is typically richer 
than the MFA phase diagram. In particular, models based 
on susceptible-infective-recovered-susceptible (SIRS) in- 
fection dynamics on lattices have been shown to exhibit 
in the PA a region in parameter space corresponding to 
stable cycles [10111112] . This region is small enough to be 
easily missed in a coarse grained numerical study but it is 
large in the admissible parameter region of an important 
class of diseases for which the susceptible turn over rate 
is much slower than the recovery rate (the weak driving 
regime). In particular, this oscillatory phase corresponds 
to realistic parameter estimations for childhood infectious 
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diseases and could be related with the phenomenon of 
non-seasonal recurrent epidemics found in different data 
records in the pre-vaccination period for this class of in- 
fections |2l9l2()j . 

Predator-prey systems are another class of population 
dynamics where non-seasonal cycling behavior has been 
reported for long 21 22 23 , and robust simple models ex- 
hibiting unforced oscillations have been much sought for 
11115116124125126127128125] . In this case also, sustained os- 
cillations have been identified in a small region of the 
phase diagram of the PA equations of a two parameter 
predator-prey model on the square lattice |15ll6j . Not 
much is known about predator-prey interaction param- 
eters , 3(X , and this small oscillatory phase may correspond 
to biologically realistic parameter values for some predator- 
prey pairs. 

In this paper we consider these two paradigms in pop- 
ulation dynamics: an epidemic model (Section [2]) and a 
predator-prey model (Section [3]). We discuss the relevance 
of the oscillatory phase in the PA phase diagrams to ex- 
plain the cycling behavior found in nature by comparing 
the behavior of the PA deterministic equations in the os- 
cillatory phase with the results of simulations on regular 
random graphs (RRGs) of different degrees for both sys- 
tems. The PA is exact on the Bette lattice, and RRGs 
are locally tree-like, so one would expect stochastic dy- 
namics running on RRGs to be well approximated by the 
PA equations, at least for large system sizes. However, we 
found that, in both cases, the oscillations predicted by 
the PA are suppressed in the simulations, even for system 
sizes as large as N = 5 x 10 7 nodes. We expect that simu- 
lations on more realistic networks, where the PA performs 
worse than for RRGs, will also fail to exhibit oscillations. 
In this sense, the oscillatory phases are an artifact of the 
PA because they disappear in finite networks. In Section^ 
we discuss the reasons why the PA fails even qualitatively 
to describe the behavior of these systems and sum up the 
conclusions of the present study. 



2 The epidemic model 

In this section, we consider the dynamics of the SIRS epi- 
demic model on a random network of homogeneous degree 
k and N nodes, a regular random graph of degree k (RRG- 
k) . A preliminary report on the results of this section can 
be found in Ref. [T3]. In the SIRS model each node can 
be occupied by an individual in susceptible (S), infected 
(/), or recovered (R) state and the following update rules 
with asynchronous update are applied. Infected individu- 
als recover at rate S, recovered individuals lose immunity 
at rate 7, and infection of a susceptible node occurs at 
infection rate A multiplied by the number of its infected 
nearest neighbors n € {0, 1, ... , k}: 
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Fig. 1. Phase diagram in the (A, 7) plane of the PA SIRS 
model for 8 = 1 and k = 2.1, 3, 4, 5 (darker shades of gray for 
smaller values of k). Region I (regions II and III) represents 
susceptible-absorbing states (active states with nonzero infec- 
tive and susceptible densities). The dashed lines correspond to 
the transcritical bifurcation curves, see Eq. (J3J); the dotted lines 
are associated with a supercritical Andronov-Hopf bifurcation. 
In region II the fixed points are asymptotically stable nodes 
or asymptotically stable foci, in region III asymptotically sta- 
ble solutions are limit cycles. The markers depicting different 
phases (I, II, III) are shown for k = 2.1. 



In the infinite population limit, with the assumptions 
of spatial homogeneity and uncorrelated pairs, the system 
is described by the deterministic equations of the standard 
or uncorrelated PA 191141: 



dP s 
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In the above equations, the variables Ps, Pi stand for the 
probability that a randomly chosen node is in state S, I, 
and the variables Psi, Psr, Pri stand for the probability 
that a randomly chosen pair of nearest neighbor nodes is 
an SI, SR, RI pair, as N —> 00. As expected, neglecting 
the pair correlations and setting the pair state probabil- 
ities equal to the product of the node state probabilities 
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these equations reduce to the classic equations of the per- 
fectly mixed or MFA SIRS model 0. 

The phase diagram of the PA SIRS deterministic model 
[Eq. (|2J|] for k — 4 was analyzed in Ref. [9]. The phase di- 
agram of the same model for several values of k in the 
range k > 2 is plotted in Fig. [TJ where we have set the 
time scale so that 5 = 1. The shades of gray of the corre- 
sponding critical curves are darker for decreasing k. In the 
following discussion we will refer to the markers depicting 
different phases (I, II, III) which are shown for k — 2.1 in 
Fig. [TJ Region I (regions II and III) represents susceptible- 
absorbing (active) states. The critical lines separating the 
absorbing and the active phases (the dashed lines) corre- 
spond to the transcritical bifurcation curves and are given 
by: 



Within the active phase, the fixed points can be asymp- 
totically stable nodes or asymptotically stable foci as in 
region II, in region III asymptotically stable solutions are 
limit cycles. The critical curves corresponding to a super- 
critical Andronov-Hopf bifurcation of the nontrivial equi- 
librium (the dotted lines) separate the active phase with 
constant densities from an active phase with stable oscil- 
latory behavior. The oscillatory phase is large for k > 2 
and it gets thinner as k increases, but it persists in the 
whole range 2 < k < 6. 

A similar phase diagram, with the Andronov-Hopf bi- 
furcation critical line bounding an oscillatory phase, was 
reported in other studies of related models on lattices 
[101111121 , where SIR dynamics with different mechanisms 
of replenishment of susceptibles is modelled at the level of 
pairs with the standard PA or another closure approxima- 
tion. All these different models exhibit an oscillatory phase 
in the regime of weak driving through slow introduction of 
new susceptible individuals (small 7 in the present case). 

The analysis of the SIRS model with the update rules 
given by Eq. (JTJ was performed by using both the PA and 
stochastic simulations on the square lattice in Ref. [T4"] . 
However, the critical line bounding an oscillatory phase 
seems to have been missed in the previous study, and for 
parameter values in the oscillatory phase no global oscilla- 
tions were found in the simulations on the square lattice. 

We have compared the behavior of the PA SIRS model 
[Eq. ©] for k = 3,4 with the results of stochastic simu- 
lations on RRGs for several system sizes. In the stochas- 
tic simulations, the system was set in a random initial 
condition with given node and pair densities, and an effi- 
cient algorithm for stochastic processes in spatially struc- 
tured systems pS] was implemented to update the states 
of the nodes according to the processes of infection, recov- 
ery and immunity waning, see Eq. ([JJ. In the nonspatial 
case, this algorithm reduces to the well-known Gillespie's 
method [32] . Unless stated otherwise, the results were ob- 
tained for RRGs with N = 10 6 nodes. In the surviving 
stochastic simulations, the system always approaches a 
well-defined nontrivial equilibrium in an over-damped or 
under-damped way, exhibiting small fluctuations around 
the average trajectory. To compare with the solutions of 



the deterministic model, for each set of parameter values 
and initial conditions, the simulations were averaged over 
10 3 realizations of a RRG. 

The results of some of these stochastic simulations on 
a RRG-4 and solutions of the PA SIRS equations [Eq. ©] 
are shown in Fig. [5] (region II) and in Fig. [3] (region III). 
The susceptible (gray lines) and the infective (black lines) 
densities are shown in Fig. [3] for two sets of parameter 
values with fixed infection rate A and decreasing rate of 
immunity waning 7. The numerical solutions of the PA 
SIRS equations are plotted in dashed lines, and the results 
of the simulations in solid lines. For parameter values well 
within region II of the phase diagram as in Fig. HJa) there 
is excellent agreement between the solutions of the PA 
SIRS model for the same initial densities and the results of 
the stochastic simulations, both for the transient behavior 
and for the steady states. This agreement deteriorates as 
7 decreases and the boundary of the oscillatory region is 
approached, see Fig.[5Jb). 

The susceptible (gray lines) and the infective (black 
lines) densities are shown in Fig. [3J for two sets of initial 
conditions for parameter values in the oscillatory region 
III. Most simulations (solid lines) die out after a short 
transient, see Fig. |3Ja), while the corresponding solutions 
of the PA SIRS deterministic model (dashed lines) con- 
verge to the stable limit cycle for all initial conditions (a 
typical set is chosen in the plot). By choosing initial con- 
ditions not far from the stable cycle predicted by the PA 
SIRS model to avoid extreme susceptible depletion dur- 
ing the transient, damped oscillations towards a nontriv- 
ial equilibrium may also be observed in region III. In Fig. 
[3Jb) a plot is shown of one of these surviving simulations 
(solid lines), together with the solution of the PA equa- 
tions (dashed lines) for the same parameter values and 
initial conditions. Thus, instead of an oscillatory phase, 
the stochastic model on a RRG-4 exhibits in region III a 
bistability phase, even for large system sizes. 

As is seen in Fig. [JJ the oscillatory phase for the PA 
SIRS deterministic model for k = 3 is large in compari- 
son with that for k = 4. It is easier in this case to avoid 
stochastic extinctions and to compare with the numerical 
simulations on a RRG-3 for a broad range of parameter 
values. Nevertheless, we have found that in the endemic 
phase of the phase diagram the stochastic dynamics of the 
SIRS on a RRG-3 is qualitatively the same as on a RRG-4. 
For large 7 it is well approximated by the PA determin- 
istic equations but no global oscillations are observed for 
small 7. 

This failure of the PA model to capture the qualita- 
tive behavior of the simulations on large RRGs has been 
investigated. Extinctions due to finite size are one of the 
reasons why the oscillatory phase is seen as an absorbing 
phase in the stochastic simulations. Indeed, as can be seen 
in Fig.[3Ja), the oscillations predicted by the PA SIRS de- 
terministic model attain very small densities of infectives 
during a significant fraction of the period in the transient 
regime. For example, the density of infectives attains val- 
ues smaller than 10~ 5 in the transient regime for the initial 
conditions of Fig. Ella) . This problem can be overcome ei- 
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Fig. 2. For 5 = 1, k = 4, comparison of the solutions of the PA deterministic model (dashed lines) with the results of stochastic 
simulations (solid lines) on a RRG-4 with N = 10 nodes for parameter values in region II. Susceptible (infective) densities are 
plotted in gray (black). Parameters: (a) 7 = 2.5, A = 2.5; (b) 7 = 0.1, A = 2.5. 
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Fig. 3. For 8 = 1, fe = 4, 7 = 0.025, A = 2.5, comparison of the solutions of the PA deterministic model (dashed lines) with the 
results of stochastic simulations (solid lines) on a RRG-4 with N = 10 6 nodes for parameter values in region III. Susceptible 
(infective) densities are plotted in gray (black). Initial conditions: (a) P s (0) ~ 0.9240, P/(0) « 0.0731, Ps/(0) w 0.0558, 
Psh(0) « 0.0024, Pri(0) w 0.0005; (b) P s (0) w 0.4889, Pj(0) w 0.0287, P S /(0) w 0.0104, P S h(0) w 0.2369, Pr/(0) w 0.0129. 



ther by setting a RRG in a random initial condition with 
node densities lying very close to the limit cycle and/or 
by increasing system size. Though at low rate of immunity 
waning the SIRS dynamics on RRGs is blurred by strong 
fluctuations and stochastic extinctions, it is quite surpris- 
ing that global oscillatory behavior has not been identified 
in stochastic simulations at all. Increasing system size up 
to N = 5 x 10 7 nodes we still find suppression of oscil- 
lations in region III and significant discrepancies between 
the transient and steady states of the PA SIRS solutions 
and the results of the simulations in region II close to the 
boundary with region III. 



3 The predator-prey model 

A small oscillatory phase has also been identified in the 
phase diagram of the PA equations of a two parameter 
predator- prey model on the square lattice [1.5116] . In this 
section we consider a generalization of this model and 
study whether the oscillations predicted by the PA equa- 
tions persist in stochastic simulations on RRGs. 

Consider a predator-prey model in which each node of 
a RRG-fc can be either empty (E) or occupied by a preda- 
tor (A) or a prey (B). Let ni,n,2 € {0,1,..., fc} denote 
the number of nearest neighbors occupied by a prey in the 
neighborhood of a node in state E and by a predator in the 
neighborhood of a node in state B, respectively. The dy- 
namics of the system at the microscopic level is governed 
by the following four processes: death of predators with 
rate d, birth of prey with rate b ■ n\ , and two competing 
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Fig. 4. Phase diagram for the predator-prey model described by Eq. (JS| for k = 4. The constraints on the independent 
parameters p and d given by Eq. (J7J) define the triangle area as parameter space for fixed k ■ p2- (a) Region I represents prey- 
absorbing states and regions fl and III represent active states with nonzero densities of predators and prey. Several critical lines 
of transcritical bifurcations between prey-absorbing and active phases are shown for k = 4 and different values of P2, see Eq. 
l|8) . (b) The active states can be asymptotically stable nodes or asymptotically stable foci as in region II or stable limit cycles as 
in region III. The phase with stable oscillatory behavior is bounded by supercritical Andronov-Hopf bifurcation curve plotted 
for k — 4 and different values of p2 for comparison. 



predator-prey interactions with rates p\ ■ ri2 and P2 ■ ni as- 
sociated with predator reproduction and consumption of 
prey. The state of the system evolves in time with asyn- 
chronous update according to the following set of local 
rules: 



(4) 



E 
B> 



b-ni 



B P2 ^ 2 E. 



The standard PA model for this dynamics on a RRG 
of coordination number k is described by: 



dP A 
dt 

dP B 
dt 

dPAB 

dt 



dPAi 
dt 



= -dP A + kp 1 P AB , 

= -k{pi +p 2 )Pab + kbP B E , 

I j. j.^p i (k-l)bP AE PBE 
= -(P1+P2 + d)P AB -' 
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PiPafs ^ „ (PiP B e +P2PAB) , 



dP BE _ (fc - l)P A i 



dt P B 

, (k-l)bP BE 



1-P A -P B 
bP BE + dP AB , 



Pb 

[p 2 (P B - Pas - 2Pbe) - PiPbe] 
1-Pa-Pb-Pae- 2P be ) ~ 



where, using the notation of the previous section, the vari- 
ables stand for the limit values of the node P A , P E and 
pair densities P AE , P A e, Pbe as N — > 00. 

Linear stability analysis of the equilibrium points of 
Eq. ([5]) allows to identify different regions in the phase 
diagram, whose parameters are the coordination number 
of a RRG and the four rate constants of the model. For 
Pi = 0, Eq. ([5]) reduces to the model studied in Ref. [To] . 
Following [15] , the number of independent parameters can 
be reduced by taking 



^ = ( \ 



kb = ( - - p - - 



(6) 



so that k{p\ + b) + d = 1, leaving two independent param- 
eters p and d with the restrictions 



1/2 <p < 1/2 , < d< 1 



(7) 



For fixed fc, the parameter space is then given by a right 
triangular prism in which the base is restricted to the al- 
lowed values of (p, d) and the height is defined by the 
parameter fc ■ p 2 > 0. Equivalently, the phase diagram 
can be plotted in the (p, d) plane for the intersections 
fc • P2 = const of the prism. 

In Fig. [4] we have analyzed the phase diagram of the 
predator-prey model described by Eq. (|5|) for fc = 4. In Fig. 
2£a) region I represents prey-absorbing states in which a 
trivial fixed point is stable, regions II and III correspond 
to active states containing nonzero predator and prey den- 
sities. For fc = 4 and fixed P2, the critical lines d c = f(p) 
separating prey- absorbing states from active states corre- 
spond to the transcritical bifurcation curves that are given 
by the positive root of the quadratic equation: 



kp2 



3 - 64 - 29rfg + 32d c p - 12p 2 
8d c 



(8) 
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Fig. 5. For k — 4, p — 0, d = 0.017, P2 = 0, comparison of the solutions of the PA deterministic model given by Eq. (J5]) (dashed 
lines) with the results of stochastic simulations (solid lines) on a RRG-4 with N = 10 6 nodes for parameter values in region III. 
Predator (prey) densities are plotted in gray (black). The three panels show time interval of equal length in the beginning, in 
the middle and in the end of the same time series. 



Several solutions of Eq. ((HJ for k = 4 and different values of 
P2 are shown in Fig. 3(a). As for the asymptotic behavior 
of the system within the active phase, Fig. [4(b) shows two 
different regions: region II correponds to asymptotically 
stable nodes or asymptotically stable foci (not separated 
in the plot) and region III represents stable limit cycles. 
The stability of the fixed point in region II is lost on the 
critical line separating regions II and III through a super- 
critical Andronov-Hopf bifurcation where a limit cycle is 
born. Again, several critical curves for k — 4 and different 
values of P2 are plotted in Fig. 1Kb). Note that the oscil- 
latory phase occuring for small values of d is also present 
for other sets of parameters, in particular for k ^ 4, but in 
general it is small when compared with the whole parame- 
ter space (see the scale in the left and in the right panels in 
Fig-H]for comparison). For fixed P2, the dependence of the 
oscillatory phase on k is similar to that of the SIRS model. 
Moreover, the oscillatory phase is robust with respect to 
several perturbations of the model (for example, with re- 
spect to the inclusion of the process describing natural 
death of prey). The analysis of a simplified model with 
P2 = and k = 4 was performed by using both the PA 
and stochastic simulations on the square lattice [TS]. For 
parameter values in the oscillatory phase, no global oscil- 
lations were found in the simulations. More surprisingly, 
we have found that the oscillations are also suppressed in 
simulations on a RRG-4, even for system sizes as large as 
N = 5 x 10 7 nodes. 

The typical results of the averaged stochastic simula- 
tions on a RRG-4 and solutions of the PA model given by 
Eq. ([5]) for the same initial conditions and parameter val- 
ues in the oscillatory region III are shown in Fig. [5) The 
figure is composed by three panels showing time intervals 
of equal length in the beginning, in the middle and in the 
end of a long time series. The predator (prey) densities 
are plotted in gray (black) . The numerical solutions of the 
PA equations (dashed lines) converge to a limit cycle but 



the results of the simulations (solid lines) exhibit slowly 
damped oscillations towards a nontrivial equilibrium. 



4 Discussion and conclusions 

The breakdown of the PA in the oscillatory phase found 
for both the epidemic and the predator-prey model is 
partially due to stochastic extinctions in the simulations. 
However, by choosing carefully the initial conditions, the 
simulations persist for long times exhibiting damped os- 
cillatory behavior towards a nontrivial equilibrium, which 
shows that there are other effects at play. 

We shall discuss these in the framework of the SIRS 
model, which is simpler, although a similar reasoning ap- 
plies to the predator-prey model. As pointed out in Sec- 
tion [2 a crucial parameter for the performance of the PA 
against stochastic simulations is the rate of immunity wan- 
ing 7. In the limit 7 — s- 0, the SIRS model coincides with 
the SIR model, in the opposite limit 7 — s- 00, with the SIS 
model, also known as the contact process [216117] . At high 
7 the process of loss of immunity by recovered individuals 
dominates over the infection that proceeds via SI pairs. 
In this regime pairs are constantly changing type. Due to 
such randomization or homogenization of states and to 
the randomness of the connections between the nodes in a 
RRG, the correlations are easily destroyed and an effective 
mixing of the population is achieved. Thus, it is expected 
that, in the limit 7> 1, the inclusion of short-range cor- 
relations in a deterministic model is sufficient for accurate 
description of the results of stochastic simulations. 

In the limit 7 <§; 1, the dynamics of the model is drasti- 
cally different. In this regime the rate of immunity waning 
is smaller than the recovery rate S. As pointed out in the 
introduction, the corresponding region in the phase dia- 
gram is related in the epidemiological context with acute 
disease spread 19. . Recovery and infection governed by SI 
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pairs dominate when 7 <C 1, and a RRG is crowded with 
recovered individuals while the overall number of infec- 
tives remains low. The density of SI pairs through which 
the infection spreads is even lower. In this regime both 
stochastic effects due to finite size and any kind of " imper- 
fections" in the structure of the RRG become important. 
Indeed, the standard PA is exact for tree-like structures 
where each node has exactly the same number of neighbors 
and there are no loops, the Bethe lattices. These infinite 
structures cannot be simulated on a computer. On the 
other hand, loops are inherent to all cZ-dimensional spatial 
structures used in computer simulations. Classic results 
of graph theory show that a particular realization of a 
RRG-A: will contain a large number of loops, of which the 
overwhelming majority are long (with respect to the av- 
erage path length), so that locally the graph is essentially 
tree-like. One would expect then the qualitative predic- 
tion of the PA to perform well on RRGs, provided they 
are large enough. It is striking how the dynamics of the 
SIRS unfolding on RRGs appears to be so different from 
that predicted by the PA deterministic equations, reveal- 
ing the subtle influence of long loops and stochastic effects 
whose role increases with decreasing k. 

We relate the observed suppression of oscillations also 
with the fact that in the PA model the oscillatory phase 
disappears for large k, where we could expect the PA to 
perform better. A similar observation of emergence and/or 
suppression of global oscillations in the qualitative and 
quantitative comparison of the results of Monte Carlo sim- 
ulations on RRGs with the predictions of the standard 
PA was reported in Refs. [55J and [53] for a spatial Rock- 
Scissors-Paper game, a system where three states cycli- 
cally dominate each other. In this model, for k > 4 the 
numerical solutions of the standard PA equations exhibit 
increasing amplitude oscillations of the state probabilities. 
Stochastic simulations of the model on RRGs show differ- 
ent qualitative behavior. For k = 3, 4 the evolution tends 
towards a limit cycle and the growing spiral trajectories 
are observed on RRGs only for k > 6. The six-node ap- 
proximation predicts a limit cycle in good quantitative 
agreement with the simulations on a RRG-3 |33I34| . 

The standard PA considered in this paper can be ap- 
plied straightforwardly for modeling the corresponding dy- 
namics on networks with heterogeneous topologies. In this 
case, the parameter k would stand for the first moment of 
the degree distribution. However, our results for random 
networks with homogeneous degree lead us to expect an 
even worse agreement with the numerical simulations on 
general heterogeneous networks. A number of recent stud- 
ies have indicated that an improvement in the description 
of simple dynamics on more realistic network structures 
might be possible by introducing explicitly the degree de- 
pendence of the node and pair probabilities in the PA 
equations |17I18I35] . 

In conclusion, we have shown that two classes of mod- 
els in population dynamics, an epidemic and a predator- 
prey model, with the coordination number k as a param- 
eter, exhibit similar phase diagrams in the PA. A distinc- 
tive feature of these diagrams is the presence of a small 



phase with stable oscillatory behavior in a weak driving 
regime of the active phase that is bounded by a curve as- 
sociated with a supercritical Andronov-Hopf bifurcation. 
This oscillatory phase changes with k and vanishes for 
large k in both models. The oscillatory phase has been 
associated with the possibility of unforced stable oscilla- 
tions in predator-prey models and non-seasonal recurrent 
epidemics, however, even for RRGs, for which the PA is 
known to perform quite well, the results of stochastic sim- 
ulations do not confirm the analytic prediction. In this 
regime the numerical results demonstrate that apart from 
finite size effects the macroscopic quasi-stationary state 
of the system is a nontrivial fixed point. The cyclic global 
behavior predicted by the PA is not observed in finite 
systems up to very large sizes. In order to capture the 
dynamics of these models in the weak driving regime, an 
analytic description based on more elaborate approxima- 
tion schemes than the PA must be considered, even when 
the underlying interaction network is a random graph. 
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